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ABSTRACT 

Certain  new  approaches  to  linear  programming  have  recently  received  considerable  publicity  be¬ 
cause  of  the  promise  of  substantial  improvements  in  efficiency  compared  to  the  simplex  method. 
This  note  briefly  discusses  several  research  directions  in  methods  for  solving  linear  programs  using 
nonlinear  problem  transformations.  In  particular,  we  describe  application  of  a  barrier  transforma¬ 
tion  to  the  dual,  and  the  development  of  sparse  least-squares  methods  based  on  the  LU  factorization 
of  the  least-squares  matrix  or  its  transpose. 


The  material  contained  in  this  report  is  based  upon  research  supported  by  the  U.S.  Department 
of  Energy  Contract  DE-AA03-76SF00326.  PA  No.  DE  AS03-76ER72018;  and  the  Office  of  Naval 
Research  Contract  N0001 1-85-K-0343. 

t  Presented  at  the  Workshop  on  Iterative  Methods  for  Mathematical  Programs,  Stanford  Univer¬ 
sity,  April  17  18,  1986. 


,•  •  -  m 


2.  A  Barrier- Function  Method  for  the  Primal  LP 


1.  Introduction 


Linear  programming  (LP)  has  been  the  subject  of  great  interest  recently  because  of  publicity 
surrounding  an  algorithm  (Kanuarkar,  1984)  that  is  not  only  polynomial  in  complexity,  but  also  is 
claimed  to  be  much  faster  than  the  simplex  method,  developed  nearly  40  years  ago  by  George  B. 
Dantzig.  Since  the  publication  of  Karmarkar’s  original  algorithm,  interest  has  continued  to  grow  in 
alternatives  to  the  simplex  method  (see,  e.g.,  Gay,  1985;  Goldfarb  and  Mehrotra,  1985;  Osborne, 
198G;  Sltanno  and  Marsten,  1985;  Todd  and  Burrell,  1985;  Tomlin,  1985;  Vanderbei,  Meketon  and 
Freedman,  1985).  In  this  connection,  Gill  et  al.  (1985)  have  shown  that  Karmarkar’s  projective 
method  is  closely  related  to  the  classical  logarithmic  barrier  function.  In  particular,  Karmarkar’s 
projective  method  is  a  special  case  of  a  projected  Newton  method  applied  to  the  logarithmic 
barrier  function.  Based  on  this  equivalence,  many  researchers  are  re-examining  the  effectiveness  of 
once-discarded  methods  that  in  effect  “nonlinearize”  a  linear  program.  The  proof  of  polynomial 
complexity  of  Karmarkar’s  original  algorithm  suggests  that  a  suitable  nonlinear  transformation 
may  somehow  overcome  the  inherently  combinatorial  nature  of  the  simplex  method. 

Because  the  main  interest  in  nonlinear  LP  techniques  arises  from  their  potential  speed,  an 
essential  part  of  the  work  reported  in  Gill  et  al.  (1985)  involved  numerical  experiments  with  an 
initial  implementation  of  the  barrier  method  on  a  set  of  small-  to  medium-scale  test  problems.  The 
results  were  encouraging  in  that  the  barrier  method  was  comparable  in  speed  to  the  simplex  method 
on  certain  problems.  Typically,  a  barrier  algorithm  converges  in  a  smaller  number  of  iterations 
than  the  simplex  method.  However,  since  each  iteration  of  the  barrier  method  requires  solution  of 
a  linear  least-squares  subproblem,  a  barrier  method  will  be  faster  than  the  simplex  method  only 
if  these  subproblcms  can  be  solved  quickly.  Thus,  the  hope  is  that  barrier-type  methods  will  be 
faster  than  the  simplex  method — even  substantially  so — for  problems  with  structure  that  allows 
fast  solution  of  the  least-squares  subproblems. 

The  remainder  of  this  note  is  organized  as  follows.  In  Section  2,  we  summarize  some  back¬ 
ground  on  the  application  of  a  logarithmic  barrier  function  method  to  a  linear  program  in  standard 
form.  A  disadvantage  of  the  primal  formulation  is  that  the  search  direction  is  restricted  to  lie  in  a 
subspacc.  In  order  to  satisfy  this  requirement,  cither  the  associated  linear  least-squares  subproblem 
must  be  solved  accurately,  or  the  vector  must  be  projected  into  the  appropriate  subspace.  As  an 
alternative,  Section  3  gives  the  dual  formulation  of  an  LP,  which  leads  to  a  purely  unconstrained 
subproblem,  thereby  allowing  approximate  methods  to  be  used  to  solve  the  least-squares  subprob- 
lern.  In  Section  4,  we  discuss  techniques  for  solving  the  least-squares  subproblem,  and  consider 
possible  uses  of  the  LU  factorization. 

2.  A  Barrier-Function  Method  for  the  Primal  LP 

This  section  is  intended  to  serve  as  general  background  on  barrier  functions  and  on  the  work 
described  in  Gill  ct  al.  (1985). 

Barrier-function  methods  treat  inequality  constraints  by  creating  a  barrier  function ,  which  is 
a  co  nhination  of  the  original  objective  function  and  a  weighted  sum  of  functions  with  a  positive 
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singularity  at,  tilt-  constraint  boundary.  (Many  harrier  functions  have  been  proposed;  we  consider 
only  the  logarithmic  barrier  function,  first  suggested  by  Frisch,  1955.)  Barrier-function  methods 
require  a  strictly  feasible  starting  point  for  each  minimization,  and  generate  a  sequence  of  strictly 
feasible  iterates.  (For  a  complete  discussion  of  barrier  methods,  see  Fiacco,  1979;  both  barricr- 
and  penalty-function  methods  are  described  in  Fiacco  and  McCormick,  1968.  Brief  overviews  are 
given  in  Fletcher,  1981,  and  Gill,  Murray  and  Wright,  1981.) 

In  particular,  consider  the  inequality-constrained  problem 

minimize  <f>(x) 

*eRn  (2.1) 

subject  to  h(x)  >  0, 

where  h  is  an  l-vector  of  nonlinear  functions.  Let  x  denote  the  solution  of  (2.1).  With  a  barrier- 
function  approach,  the  inequality  constraints  of  (2.1)  are  “removed”  by  including  transformed 
versions  in  a  modified  objective  function,  as  follows.  Given  a  barrier  parameter  p  (p  >  0),  let  x(p) 
denote  the  unconstrained  minimum  of 


B(x,p)  =  <f>(x)  -P^,  ln(M*))- 


Under  mild  conditions,  it  can  be  shown  that 

lim  x(p)  =  x  . 
o 

Many  results  have  been  proved  concerning  the  asymptotic  properties  of  the  sequence  (a:(^c)}  is 
p  — ►  0  (see,  e.g.,  Mifflin,  1972a,  b;  Jittorntrum,  1978;  Jittorntrum  and  Osborne,  1978). 

Now  consider  applying  a  barrier-function  method  to  the  following  linear  program  in  standarH 
form  (which  will  also  be  called  the  primal  LP): 

minimize  cTx  (2.2a) 

subject  to  Ax  =  b  (2.2 b) 

x  >  0,  (2.2c) 

where  A  is  an  m  x  n  matrix  with  m  <  n. 

Since  the  barrier  transformation  may  be  applied  only  to  the  inequality  constraints  (2.2c), 
the  subproblem  associated  with  the  primal  LP  treats  the  linear  equality  constraints  directly  and 
transforms  only  the  bounds: 

n 

minimize  F(x)  =  cTx  -  p)  lnx< 

»e*"  (2.3) 

subject  to  Ax  =  6. 

A  standard  approach  to  solving  a  linearly  constrained  problem  of  the  form  (2.3)  is  to  use  a 
feasible-point  descent  method  (see,  e.g.,  Gill,  Murray  and  Wright,  1981).  The  current  iterate  x 
always  satisfies  Ax  =  b,  and  the  next  iterate  x  is  defined  as 

x  =  x  +  ap,  (2.4) 
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when1  p  is  an  n-vector  (the  search  direction),  and  a  is  a  positive  scalar  (the  stcplcngth).  The 
computation  of  p  and  a  must  ensure  that  Ax  =  b  and  F(x)  <  F(x). 

The  Newton  search  direction  is  defined  as  the  step  to  the  minimum  of  the  quadratic  approx¬ 
imation  to  F(x)  derived  from  the  local  Taylor  series,  subject  to  retaining  feasibility.  Thus,  the 
Newton  search  direction  is  the  solution  of  the  following  quadratic  program: 

minimize  gTp  +  \pTHp  (2.5a) 

Pfc»“ 

subject  t«  Ap  —  0,  (2.56) 

where  g  =  VF(x)  and  H  =  V2F( x).  The  optimality  conditions  for  (2.5)  imply  that 

g+  Hp-  Atx  (2.6) 

for  a  vector  n  of  Lagrange  multipliers.  Thus,  p  and  ir  satisfy  the  partitioned  linear  system: 

c  :)(-;)=(;)• 

if  #  is  nonsingular,  one  method  of  solving  (2.7)  is  to  solve  the  equations 

AHlATx  =  AHlg ,  (2.8a) 

Hp  —  ATir  —  g.  (2.86) 

Let  pP  (the  Newton  direction  for  the  primal  LP)  denote  the  solution  of  subproblem  (2.5)  when 
F  is  the  barrier  function  in  (2.3),  with  associated  Lagrange  multiplier  *>.  The  derivatives  of  F 
are: 

g(x)  =  c-pX~le  and  H(x)  =  /iX~3,  (2.9  a) 

where 

X  =  diag  (*,),  j  =  l,...,n,  (2.96) 

and  e  —  (1, 1, . . . ,  1)T.  Since  H(x)  is  positive  definite  when  x  >  0,  pP  is  finite  and  unique,  and  is 
a  descent  direction  for  F(x),  i.e.,  (c  -  p.X~le)TpP  <  0.  Substituting  from  (2.9a)  in  (2.8 a),  we  see 
that  irP  satisfies 

AX*At kp  =  AX{Xc  -  pe).  (2.10) 

Recall  from  standard  linear  algebra  (e.g.,  Stewart,  1973)  that,  for  any  matrix  C,  the  solution 
of  the  least-squares  problem 

minimize  \\f  -  Cx ||2  (2.11) 

satisfies  the  normal  equations 

CTCx  =  CTf ,  (2.12) 

which  are  always  compatible.  In  referring  to  (2.11)  and  (2.12),  we  shall  call  C  the  least-squares 
matrix  and  CTC  the  normal-equation  matrix. 
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Comparing  (2.10)  and  (2.12)  and  taking  C  —  XAT  in  (2.12),  it  follows  that  7r,.  solves  the 


least-squares  problem 

minimize  \\Xc  —  fie  -  XATn jj*. 

(2.13) 

The  vector  pP  is  then  defined  by 

pr  =  -(1  /p)XrP, 

(2.14) 

where  rP  is  the  optimal  residual  of  (2.13): 

tp  —  Xc  —  fie  —  XAttcp. 

(2.15) 

Note  that  if  m  variables  are  nonzero  at  the  solution  of  (2.2)  (i.e.,  if  the  primal  is  nondegenerate), 
then  XAt  will  retain  full  rank  even  as  the  iterates  approach  the  solution.  (This  contrasts  with 
the  structural  ill-conditioning  of  the  Hessian  matrices  of  the  barrier  function  when  fewer  than  n 
constraints  are  active  at  the  solution;  see,  e.g.,  Murray,  1971.)  However,  a  well  known  feature 
of  practical  linear  programs  is  that  they  usually  arc  degenerate;  i.e.,  fewer  than  m  variables  are 
nonzero  at  the  solution.  In  this  case,  the  matrix  XATis  rank-deficient  at  the  solution,  and  becomes 
increasingly  ill-conditioned  as  //.  — *  0.  Therefore,  any  method  used  in  practice  to  solve  (2.13)  must 
be  able  to  cope  with  singularity  and  extreme  ill-conditioning.  (Numerical  solution  of  the  least- 
squares  subproblcms  will  be  discussed  in  Section  4.) 

3.  Solving  the  Dual  Problem 

With  the  approach  described  in  Section  2,  the  search  direction  is  required  to  lie  in  a  particular 
subspace  in  order  for  the  iterates  to  retain  feasibility  with  respect  to  the  equality  constraints  (2.2 b). 
In  this  section,  we  describe  an  alternative  based  on  the  dual  linear  program  associated  with  (2.2): 

minimize  bTy 

(3.1) 

subject  to  ATy  >  —c. 

The  optimal  vector  y*  for  (3.1)  is  the  negative  of  the  Lagrange  multiplier  vector  r  for  the  equality 
constraints  of  (2.2).  (Use  of  the  dual  LP  in  this  context  has  been  suggested  independently  by  other 
researchers,  for  example  at  the  recent  ONR-sponsored  Workshop  on  New  Directions  in  Mathemat¬ 
ical  Programming,  Monterey,  California,  February  20-  21,  1986.  It  is  also  closely  related  to  the 
work  of  Eriksson,  1981,  1985.) 

Since  all  the  constraints  of  (3.1)  are  inequalities,  applying  the  barrier  transformation  to  this 
problem  gives  a  purely  unconstrained  subproblem: 

n 

minimize  F(y)  =  bTy  -  ln(ajy  +  c^).  (3.2) 

V  i= i 

Let  y(p)  denote  the  unconstrained  minimum  of  (3.2).  Defining 

D  =  diag(dj)  and  X  =  diag(zy)  =  fiD~ *, 


4 


3.  Solving  the  Dual  Problem 


where  dj  =  ajy  +  Cj  and  x j  =  fi/dj,  we  have 

VF  =  b  -  Ax  and  VJF  =  -AX'A*. 

P 

The  gradient  of  F  must  vanish  at  y(/t),  which  implies  that 

n 

afy(p)  +  Cj  °J  -  6’ 


and  hence 


M  * 

um  —  =  x.  . 

M-0  dj  1 


(3.3a) 


(3.36) 


In  order  to  define  a  Newton  method  for  solving  (3.2),  we  assume  that  a  point  y  is  available  for 
which  ATy  +  c  >  0  (so  that  Vs  F  is  positive  definite).  Using  the  above  derivatives,  pD  (the  Newton 
search  direction  for  the  dual)  then  satisfies  the  system 

AX*ATPn  =  p(Ax  -  6).  (3.4) 

In  contrast  to  the  subproblem  (2.5)  associated  with  the  primal  problem,  pD  need  not  lie  in  any 
particular  subspacc,  and  hence  (3.4)  may  be  solved  only  approximately — for  example,  using  a  few 
iterations  of  a  conjugate-gradient  method. 

Note  that  if  the  right-hand  side  of  (3.4)  can  be  expressed  in  the  form  AXv  for  some  vector 
v,  (3.4)  will  have  the  form  (2.12)  of  the  normal  equations  for  a  least-squares  problem.  (It  is 
advantageous  numerically  to  solve  a  least-squares  problem  rather  than  equations  of  the  form  (3.4). 
However,  note  that  (3.4)  always  has  a  bounded  solution  if  a  primal-feasible  point  exists.)  We 
therefore  seek  v  such  that 

AXv  =  Ax  —  6,  (3.5) 

whereupon  pD  becomes  the  solution  of  the  least-squares  problem 

minimize  ||/jv  -  Jfi4Tp||J,  (3.6) 

which  may  be  solved  using  the  techniques  to  be  described  in  Section  4.  Note  that  x  >  0,  whereas 
A(x  -  Xv)  =  6.  As  convergence  occurs,  v  — »  0  and  x  — *  x* 

It  might  appear  that  the  elements  of  the  least-squares  matrix  would  become  unbounded  with 
this  approach,  in  contrast  to  the  primal.  However,  the  relationship  (3.3)  shows  that  dj  — »  0  only 
as  p  — ►  0,  and  thus  the  limiting  matrices  X  and  XAT  should  remain  bounded. 

An  alternative  to  the  exact  dual  formulation  (3.4)  is  to  retain  the  same  “Hessian”  for  several 
iterations.  In  Newton-type  methods  for  nonlinear  optimization,  it  is  common  to  define  the  search 
direction  with  a  positive-definite  approximation  to  the  Hessian  matrix.  (For  example,  the  exact 
Hessian  may  be  indefinite  or  expensive  to  compute.)  In  the  present  context,  any  positive-definite 
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X  may  be  used  in  (3.4)  (3.6)  to  generate  a  descent  direction.  It  is  essential  to  retain  the  definition 
x  —  jiD  xc  to  obtain  the  correct  gradient  b  —  Ax,  but  X  need  hot  be  defined  as  pD  *. 

An  immediate  implication  is  that  same  matrix  X  could  be  used  for  perhaps  several  iterations 
of  the  Newton-type  method.  (Convergence  results  for  methods  for  this  type  are  given  in,  e.g., 
Dennis,  1970,  and  Ortega  and  Rheinboldt,  1970.)  Since  a  sparse  factorization  of  XylTis  needed  to 
compute  v  and  p  in  (3.5)  (3.6),  considerable  economies  can  be  made  by  retaining  the  factorization 
from  the  previous  iteration. 

4.  Solving  the  Least-Squares  Subproblem 

It  should  be  emphasized  that  the  least-squares  subproblems  associated  with  the  primal  and  dual 
linear  programs  ((2.13)  and  (3.6))  are  very  similar,  and  it  is  crucial  that  they  be  solved  efficiently. 
In  this  section  we  review  techniques  that  have  been  suggested  or  used,  and  then  describe  some  new 
approaches  based  on  the  LU  factorization. 

4.1.  Background 

The  major  strategics  for  solving  sparse  least-squares  problems  include: 

•  direct  methods,  based  on  a  sparse  QR  factorization  of  C  (e.g.,  George  and  Heath,  1980)  or  on 
the  Cholesky  factorization  of  the  normal-equation  matrix  CTC  (e.g.,  George  and  Liu,  1981); 

•  updating  methods  that  factorize  most  of  C  or  CTC  and  deal  with  a  few  omitted  rows  or 
columns  by  partitioning  (e.g.,  Heath,  1984); 

•  purely  iterative  methods  that  perform  no  factorization,  such  as  a  straightforward  conjugate- 
gradient  method; 

•  hybrid  methods  that  combine  factorization  and  iterative  techniques.  For  example,  in  Gill  et 
a 1.  (1985),  the  subproblem  (2.13)  is  solved  with  a  preconditioned  version  of  the  stabilized 
conjugate- gradient  method  LSQR  (Paige  and  Saunders,  1982),  where  the  preconditioner  is 
the  Cholesky  factor  of  a  sparse  matrix  related  to  AX2AT.  (The  use  of  a  preconditioned 
conjugate-gradient  method  in  this  context  is  also  discussed  by  Gay,  1985.)  For  early  work  on 
preconditioners,  see  Axelsson  (1974)  and  Concus,  Golub  and  O’Leary  (1976).  Some  techniques 
for  computing  an  approximate  (“incomplete”)  Cholesky  factor  based  on  sparsity  considerations 
arc  given  in,  e.g.,  Meijerink  and  van  der  Vorst  (1977),  Mantcuffcl  (1980)  and  Munksgaard 
(1980). 

The  effectiveness  of  these  options  varies  with  the  problem. 

When  using  an  iterative  scheme,  it  is  desirable  to  be  able  to  use  only  an  approximate  solu¬ 
tion  (obtained,  say,  from  a  truncated  conjugate-gradient  method;  see,  e.g.,  Dembo,  Eisenstat  and 
Stcihaug,  1982).  In  this  context,  a  disadvantage  of  the  primal  formulation  is  that  an  approximate 
solution  of  (2.13)  will  not  satisfy  the  equality  constraints  (2.56).  One  means  of  overcoming  this 
difficulty  is  to  project  the  search  direction  into  the  null  space  of  A,  for  example  using  a  sparse 
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representation  of  a  basis  for  the  null  space.  Such  a  matrix  may  be  constructed  from  an  LU  fac¬ 
torization  of  A  (which  would  need  to  be  computed  only  once),  or  alternatively  it  may  be  obtained 
from  an  LU  factorization  of  XAT. 

For  the  dual  formulation,  a  truncated  conjugate-gradient  method  is  most  promising  because 
an  approximate  solution  of  (3.4)  can  be  shown  to  be  a  descent  direction,  as  required.  However, 
the  speed  of  convergence  of  the  conjugate-gradient  method  is  crucial  to  practical  success.  We  now 
consider  methods  for  preconditioning  the  conjugate-gradient  method  using  an  LU  factorization. 


4.2.  The  LU  factorization 

Given  a  general  sparse  matrix  C  with  n  rows  and  m  columns,  we  write 

C  =  LU,  (4.1) 

where  L  is  square  and  nonsingular,  the  dimensions  of  U  are  those  of  C,  and  L  and  U  are  (nominally) 
unit-lower-triangular  and  upper-triangular  respectively.  (In  general,  the  rows  and  columns  of  C  will 
have  been  permuted  during  the  factorization,  but  we  shall  ignore  the  permutations  for  simplicity 
of  notation.)  If  a  suitable  form  of  threshold  pivoting  (see  Section  4.4)  is  used  in  computing  L  and 
U,  L  tends  to  remain  well-conditioned  throughout,  so  that  the  condition  of  C  is  almost  always 
reflected  in  U. 

A  set  of  procedures  (LUSOL;  sec  Gill  et  ai.,  1986)  has  recently  been  developed  for  computing 
and  updating  the  LU  factors  of  a  general  sparse  matrix;  the  techniques  in  LUSOL  are  related 
to  the  work  of  Reid  (1976,  1982)  on  computing  a  sparse  LU  factorization  and  performing  sparse 
Bartels-Golub  updates  (Bartels  and  Golub,  1969).  LUSOL  provides  a  means  for  implementing 
several  methods  to  be  described  below. 

In  the  problems  of  interest,  C  is  of  the  form  X AT  and  n  >  m.  Hence  (4.1)  has  the  form 

C  =  XAT=LU  =  LU,  Z  =  U=i^j,  (4.2) 

where  L  denotes  the  first  m  columns  of  L  and  U  is  an  m  x  m  upper-triangular  matrix.  (The  last 
( n  -  m)  columns  of  L  will  simply  be  those  of  the  identity.)  Again  L  should  be  well-conditioned, 
while  U  reflects  the  condition  of  C.  For  the  present  we  assume  that  C  has  full  rank,  so  that  0  is 
nonsingular. 

Note  that  (4.2)  can  be  used  to  find  the  vector  v  in  (3.5)  needed  to  formulate  the  dual  as  the 
least-squares  subproblcm  (3.6).  It  is  sufficient  to  solve  the  two  systems 


UTu  =  Ax  -  b 


4.3.  LU  preconditioning 

Following  Peters  and  Wilkinson  (1970)  and  Bjorck  (1976),  we  note  that  the  least-squares  problem 

minimize  ||/  -  Cp\\l  (4.3) 
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may  be  solved  using  a  factorization  of  the  form  (4.2),  by  solving 

(4.4a) 
(4.41.) 


minimize  ||/  —  £<7||f, 
? 

Up  =  q. 


Because  L  is  well-conditioned,  an  iterative  method  may  converge  more  rapidly  on  problem  (4.4a) 
than  on  the  original  problem  (4.3)  (unless  C  itself  is  well-conditioned  or  has  clustered  singular 
values).  Some  experiments  with  this  approach  are  described  in  Saunders  (1979);  for  related  work, 
see  Delves  and  Barrodale  (1979). 

In  the  present  context,  greater  efficiency  is  likely  to  result  from  further  preconditioning,  par¬ 
ticularly  as  the  solution  is  approached.  If  L  is  the  first  m  rows  of  L: 


it  follows  that  p  may  be  determined  by  solving 

(4.5a) 
(4.56) 


minimize  ||/  -  )  »lll. 

Lq  =  s ,  Up  ~  q. 


The  nature  of  the  linear  programming  problem  and  the  pivoting  strategy  used  in  LUSOL  arc  such 
that  ML~l  — *  0  as  y  — ♦  y*  Hence,  the  number  of  iterations  required  by  a  conjugate-gradient 
algorithm  should  become  small  as  the  solution  is  approached. 

In  fact,  the  factorization  XAT  ~  LU  obtained  with  the  LUSOL  pivoting  strategy  serves  to 
partition  X  and  A  in  the  form 


and  (4.5)  is  equivalent  to  solving  the  original  problem  (4.3)  with  the  preconditioner  XBDT  =  LU. 
At  all  stages,  the  partitioned  factorization  pinpoints  a  nonsingular  B  that  corresponds  to  the  usual 
basis  matrix  in  the  simplex  method,  and  ultimately  it  will  correspond  to  an  optimal  basis. 

In  general,  it  may  not  be  necessary  to  recompute  the  LU  factorization  at  every  iteration.  If 
the  current  factors  are  Xi,AT  =  Li,Uk,  the  factors  of  Xk+\AT  could  be  taken  as  =  DLk  and 
Us, +i  =  Uk,  where  D  =  Xh+1Xj^1,  as  long  as  D  is  reasonably  well-conditioned.  In  some  cases 
it  may  be  efficient  to  update  the  LU  factors  to  account  for  those  elements  of  D  that  introduce 
ill-conditioning,  using  some  of  the  update  procedures  available  in  LUSOL. 

Another  approach  is  to  compute  the  LU  factorization  of  CT,  so  that 


Ct=AX  =  LU, 


(4-6) 
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where  L  is  nousingular  anil  U  is  ail  m  x  «  trapezoid.  In  this  case,  LUSOL  will  need  to  l»c  enhanced 
to  use  a  threshold  form  of  complete  pivoting,  in  order  to  ensure  that  the  iirst  m  columns  of  U  point 
to  a  suitable  basis  matrix  D.  If  a  matrix  C  has  full  column  rank,  a  reliable  factorization  C  =  LU 
can  almost  always  be  obtained  using  Gaussian  elimination  with  threshold  partial  pivoting.  By  this 
we  mean  that  when  a  potential  pivot  is  selected  from  the  remaining  rows  and  columns,  it  must  be 
reasonably  large  relative  to  other  elements  in  the  same  column,  but  it  need  not  be  compared  with 
elements  in  other  columns.  Virtually  all  existing  LU  software  uses  this  strategy. 

When  computing  the  factorization  (4.6),  the  column  rank  will  be  drastically  deficient,  and  the 
triangular  part  of  U  will  reflect  the  true  rank  only  if  pivots  are  chosen  to  be  reasonably  large  relative 
to  all  remaining  elements;  i.e.,  only  if  threshold  complete  pivoting  is  implemented.  Such  a  strategy 
requires  keeping  track  of  the  largest  remaining  element  at  each  stage  without  excessive  overhead. 
Any  such  enhancement  to  LUSOL  would  carry  a  benefit  whether  C  or  CT  is  being  factorized,  since 
Vzl^is  expected  to  be  rank-deficient. 
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